On the chemical composition of cosmic rays of 
highest energy 



Grzegorz Wilk 

The Andrzej Soltan Institute for Nuclear Studies, Hoza 69, 00681, Warsaw, Poland 
E-mail: wilk@fuw.edu.pl 

Zbigniew Wlodarczyk 

Institute of Physics, Jan Kochanowski University, Swiijtokrzyska 15, 25-406 Kielce, 
Poland 

E-mail: zbigniew. wlodarczykOujk. kielce .pi 

Abstract. We present arguments aiming at reconciling apparently contradictory 
results concerning the chemical composition of cosmic rays of highest energy, coming 
recently from the Auger and HiRes collaborations. In particular, we argue that the 
energy dependence of the mean value and root mean square fluctuation of shower 
maxima distributions observed by the Auger experiment are not necessarily caused 
by the change of nuclear composition of primary cosmic rays. They could also be 
caused by the change of distribution of the first interaction point in the cascade. A 
new observable, in which this influence is strongly suppressed, is proposed and tested. 
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Recently two leading cosmic ray (CR) experiments, the Pierre Auger Collaboration 
(Auger) |1] and The High Resolution Fly's Eye Collaboration (HiRes) |2] published 
their most recent data on the depth of maxima of extensive air showers above 10^^ 
eV. Two apparently contradictory conclusions were presented. Whereas the Auger 
collaboration cautiously concluded that their data indicated a gradual increase of the 
average mass of incoming CR with energy, HiRes stated that their data were consistent 
with a predominantly protonic composition of cosmic rays. These results started a vivid 
discussion [3] . In this note, we propose a possible reconciliation of both results with an 
indication that, perhaps, there is no need to introduce a heavy (i.e., iron) component 
in the CR chemical composition. Namely, we indicate that these results could also be 
due to the influence of the distribution of the first interaction point in the cascade, at 
least partially. We therefore propose and test a new observable, cf. Eq. ([9]), in which 
this influence is strongly suppressed. 

With increasing energy the Auger data show almost monotonic changes in the 
chemical composition of CR changing from proton to iron for two types of observables 
considered: the mean depth of the maximum of the longitudinal development of air 
showers, {Xmax), and the shower-to-shower fluctuations, the root mean square (rms) 
a{Xmax)i see Fig. [T10. For {Xmax), a such dependence can be interpreted by allowing 
for the presence of two components in CR: iron, with relative abundance a, and protons, 
with relative abundance 1 — a: 

{Xmax) = (1 ~ 0!){Xmax)p + Ct{Xmax) Fe, (l) 

(where {Xmax)p and {Xmax)Fe denote the mean depth of shower maxima for the pure 
p and Fe CR's, respectively). However, the same reasoning applied to cr {{Xmax)) lead 
to nonmonotonic dependence on a in this case (seen as nonlinear spacing between lines 
corresponding to different values of a, cf. Fig. [T]b), 

= (1 - a)orJ + aal^ + a{l - a) {{Xmax)p - {X^max)Fef ■ (2) 

This has a maximum at 
1 



"=2 



2 • (3) 

(^{XjYiax)p {X^dx) Fe) 

This is seen in Fig. [T]d, where one observes that adding iron to protons results first (for 
small a) in increased fluctuations, and only for a quite large admixture of iron (large 
a) do they decrease towards the pure iron line. For this reason, as seen in Fig. |2l 

I In all figures presented here we use both experimental data and model predictions for pure proton and 
iron primaries following [Hd]. Using information on their values of {X)p, {X)pe, Up and ape^ one can 
deduce from Eqs.(IT]) and (21), in an univocal way, the corresponding values of observables of interest for 
a given value of the parameter a; in particular, the energy dependencies of {X^ax) and a {Xmax)- The 
differences {Xmax) — c {Xmax) are evaluated directly from the definition using the energy dependencies 
of {Xmax) and <J {Xmax) as given by the models. Notice that Auger compares their data with pure 
simulations, whereas HiRes quotes data including all detector effects and compares them to the models 
after the detector simulation. It means then that, unfortunately, both approaches cannot be compared 
directly. 
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for such a simple parametrization, experimental data with similar energy behavior lead 
to quite different chemical compositions, ranging from proton dominated for {X^ax) to 
iron dominated for a {Xmax)- 

Among possible reasons of such a discrepancy, we shall concentrate on problems 
connected with the development of the cascade, in particular on the significance of the 
depth of the first interaction. It is known (cf. |5l|6l[7]) that most charged particles in the 
shower are electrons and positrons with energies near the critical energy e originated 
from the electromagnetic subshowers initiated by photons from 7r° decay. The mean 
depth of the maximum of an electromagnetic shower initiated by a photon with energy 
is 



Here Xq ^ 37 g/cm^ is the radiation length and e = 81 MeV in the air. A nuclear- 
initiated shower consists of a hadronic core which feeds the electromagnetic component 
mainly through the production of vr". Therefore, in general, for an incident nucleus of 
mass A (including protons with A = 1) and total energy the depth of the shower 
maximum is given by 



where (Xi) is the mean depth of the interaction with maximal energy deposition into 
the shower (known also as the depth of the first interaction), K denotes inelasticity and 
(n) is related to the multiplicity of secondaries produced in the high-energy hadronic 
interactions in the cascade. When the composition changes with energy, (A) depends on 
energy and (Xmax) changes accordingly. For primary nuclei with mass number A and 
energy E, the shower is, within a good approximation, simply equivalent to a bundle 
of A nucleons with energies E/A each. In the case of primary protons in the hadronic 
cascade, there is a hierarchy of energies of secondary particles in each interaction and 
a similar (approximately geometrical) hierarchy of interaction energies in the cascade. 
In this case (n) has to be understood as some kind of effective multiplicity without a 
general straightforward definition. In addition, the inelasticity K can also change with 
energy p]. 

Now, the probability of observing the first interaction in a shower at a depth greater 
than X is 



where A denotes the interaction length (and is therefore connected to the cross section, 
in our case Xp-air = 24160/(Tp_air[g/cm^] for cross-section given in [mb]). It is tempting 
to use directly the exponential distributions of showers with large X^ax to calculate 
Xi (and the proton-air cross section). However, this can be done only in the case 
of a perfect correlation between Xmax and Xi. The fluctuations existing in shower 
development modify such a relation, leading to 




(4) 



(Xmax) = {XZx [(E/A) (K/in))]) + (Xi) 



(5) 




(6) 




(7) 
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Figure 1. (Color online) {X,nax} and a {X^ax) from [I] compared with the QGSJETII 
model [U using two components, protons and iron, cf. Eqs. (P) and respectively, 
with a denoting the relative abundance of iron nuclei. 



where A = kX, and k accounts for the way the energy dissipation takes place in the 
early stages of shower evolution; it is particularly sensitive to the mean inelasticity 
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Figure 2. (Color online) The energy dependence of relative abundance of iron in 
CR as extracted from (Xmax) and a (Xmax) shown in Fig. [1] (in the frame of the 
QGSJETII model gj). 



and its fluctuations. The factor k depends mainly on the way the energy dissipation 
takes place in the early stages of shower evolution and is particularly sensitive to the 
mean inelasticity and its fluctuations. Small fluctuations in multiplicity and K result 
in smaller k 

In our case, in the absence of internal fluctuations, all showers would develop 
between the first interaction point and the maximum in the same amount of matter, 
AX = Xmax — Xi, showing perfect correlation between the X^ax and Xi. This means 
that their distributions should in this case have exactly the same shape, but shifted by 
AX, i.e., the slope of the Xmax distribution. A, should be equal to the mean interaction 
length A. This relation will, however, be affected by some inevitable intrinsic fluctuations 
in shower development taking place after the first interactioii. In this case, because of 
fluctuations in AX, the correlation between X^ax and Xi will diminish. Roughly, one 
can write that 

a{X^ax) = cT{X,)+^[a{AX)], (8) 

where a (Xi) oc (Xi) and the function ^(cr) describes the influence of shower fluctuations 
after the first (main) interaction point. Notice that for the probability distribution 

§ Assuming similar fluctuations in multiplicity and inelasticity, a model predicting large (n) of 
secondary particles leads to a smaller overall fluctuations of the cumulative shower profile of secondaries, 
i.e., to smaller factor fc [5l[7]. 

II A comprehensive studies of the shower to shower fluctuations by means of Monte Carlo simulations 
can be found in recent works [9l [10] . 



On the chemical composition of cosmic rays of highest energy 



6 



given by Eq.dH]) the fluctuation in Xi is a (Xi) = ^Var (Xi) = (Xi) = A. However, 
in the case when Xi is interpreted as the main interaction point (in which the energy 
deposition to the shower is maximal) one would obtain a gamma distribution (instead 
of Eq. ([6])) for which a (Xi) = ^yVar{Xlj = (Xi)/-^/^, where k depends on the mean 
inelasticity, (K), and determines in which of the successive interactions of a projectile 
the energy deposition to the shower is maximal^. To summarize, because of Eq. (|5]), 
where {Xmax) = {^^x) + (^1)5 propose to use the following observable in which the 
influence of fluctuations of the first interaction is strongly suppressed, 

(X^ax) - o (X_.) = {XZx [{E/A) {K/{n))]) " iH^X)]. (9) 

To test this observable, we first plot its energy behavior in Figs. [3] and H] for the, 
respectively. Auger P and HiRes [2] data. Notice that now the HiRes data, where the 
distribution of Xmax was truncated at 2a (ut denotes truncated fluctuations), show 
similar behavior as the Auger data. Notice also (cf. Fig. [3]) that {Xmax) — o" {Xmax) 
given by Eq. ([9]) still depends on models of multiparticle production and is sensitive 
to the chemical composition of CR {p and Fe initiated showers are markedly different). 
Finally, Fig. 3 also tells us that the chemical composition cannot be the origin of the 
observation by Auger [I] that {Xmax) rises too slowly with energy and approximates the 
expectation for primary Fe nuclei. In fact, experimental data seem rather to be in fair 
agreement with the hypothesis of a proton dominant composition of the primary CR 
flux (assuming, of course, that the reference models used are roughly correct). Within 
the toy model of primary composition used before (with only two components: iron 
nuclei with relative abundance a and protons with abundance 1 — a, cf. Eqs. ([T]) and 
([2])) we can again evaluate a as given by the Auger experiment but this time from 
{Xmax) — o" {Xmax)- The result is shown in Fig. [51 For reference model QGSJETII, the 
abundance of iron is roughly independent of energy {a ~ 0.05 -^ 0.1 ) and even for the 
model EPOS v. 1.99 [TTj leading to the maximal abundance of iron, a increases with 
energy rather slowly (remaining in the interval a ~ 0.15 -i- 0.3 ). 

To summarize this part, we learn from Fig. [3] that the main contribution to the 
energy dependence of {Xmax) and cr {Xmax) observed by Auger [1] comes from (Xi). 
This, however, can be affected by two factors: the cross section ainei and inelasticity 
K (in fact, not only by its mean value {K), but also by its distribution). Roughly, 
(Xi) = A ■ K, where k determines in which of the successive interactions of a projectile 
the energy deposition to the shower is maximal. For a uniform inelasticity distribution 
(in the maximal possible interval for a given {K)), one has k ~ 1 + 1.85(0.75 — {K)). As 
shown recently in [12], if gluon saturation occurs in the nuclear surface region then ap-air 
aX E > 10^^ eV increases more rapidly with incident energy than is usually estimated. 
Although in [HI [13] we have argued for an overall decrease with energy of the inelasticity 
K, its increase at energies E ^ 10^^ eV is by no means excluded [jj. Both possibilities 

% Numerically for (K) = 0.7 one has a {Xi) ^ 0.96(Xi). 

+ Recently the role of inelasticity in high energy CR was discussed in [l^ using the percolation theory 
approach. 
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Figure 3. (Color online) {X,nax) ~ c {^max) as deduced from the Auger data and 
compared to different models [1 for showers initiated by protons and iron. Note that 
experimental data prefer a proton composition. 
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Figure 4. (Color online) The same as in Fig. |3l but for the HiRes experimental data 
(in this case data are truncated at 2a) [2]. 



require an abrupt onset of some "new physics" beyond the standard model, which would 
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Figure 5. (Color online) The energy dependence of the relative abundance of iron in 
CR as extracted from {X„iax) — {Xmax) given by the Auger experiment and shown 
in Fig. [31 



be difficult to accept. It is worth mentioning as an example the elongation rate, which 
in the case of Eq. (|5]) is given by [5l [15] 



D 
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d{X,) 
dlogE' 



(10) 



As reported in [16], Auger observes the apparently abrupt change in Diq at energy 
^ 2 ■ 10^*^ eV, which could signal some new physics. If we denote -D*o = 
d [{Xmax) — 0" (Xmax)] /dlogE, then one expects here something of the order of D^q — 
DiQ = —da{Xjnax) /dlogE = 3.5 g-cm~^, depending on the increase in cross section 
adopted (chosen from existing models and predictions). However, the Auger data above 
10^* provide the value 22 g-cm~^. It can be shown that such a large value leads to a 
very strong energy dependence of the cross section, dap^air/ dlogE = OASap^air, which 
seems at the moment to be very unrealistic and contradicts even the scenario of gluon 
saturation on the nuclear surface recently proposed in [12]. A more detailed discussion 
of this problem is outside the scope of this note. On the other hand, the center of mass 
collision energies of the order of few hundred TeV observed here are well beyond those 
to be studied in the foreseeable future at LHC. This means that CR are, most probably, 
the only future source of information on the properties of interactions at these energies 
and surprises should not be ruled out a priori. 

As mentioned above the Hires data (which are truncated at 2a) [2] show similar 
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Figure 6. (Color online) Illustration of the effect of statistics as given by Eq. (fT4|) with 
= E{X) assumed to be the same as the mean value given by the QGSJETOl 
model [17], see text for details. 



energy dependence for {X^ax) — o" (Xmax) as the Auger data, c.f. Fig. HI This indicates 
that, in both cases, the crucial factors are the tails of the Xmax distributions. For a 
small sample, the values of X^ax near the maximum of the distribution are preferred, 
and the estimated mean value (Xmax) differs from the expected value, E{Xmax)- To 
investigate whether the effect observed by Auger could be connected to small statistics, 
notice that, because the distribution of distances of the first interaction is exp(— X/A), 
when calculating (X) = -^i encounters S = J2f -^i which has a gamma 

distribution. 

For small samples one in fact observes the most probable values, which for a gamma 
distribution is equal to = X{N - 1). This means that we can expect that, for a 

sample consisting with N elements one has 

1 ^ 1 / 1 \ 

i ^ ' 

whereas the true expectation values for the exponential distribution (obtained for 
N — > oo) is E{X) = A. Therefore, for a sample of N elements the estimator {X)^ 
is biased by a value of the order of 

EiX) - (X)(^) ^ A (13) 
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In Fig. iniwe compare 

(X)(oo) - (X)(jv) ^ {X)(N) 1 . . 

W(oo) W(oo) iV ^ ^ 

with data of Auger, here (X)(oo) = -E'(X) denotes the value given by the model. Because 
the reference model (here QGSJETOl) does not exactly describe the experimental data 
(even in the lower energy region, i.e. for large values of N) we use here the simple 
formula a + b/N to describe the dependence on N. To further illustrate the significance 
of tails of Xmax distributions, we examine the truncated Xmax distribution in the interval 
(0,Xc„t), cf.. Fig. El where (Xmax), (Xmax) and (Xmax) - (j{Xmax), evaluated for 
different Xcut, are presented. Notice that (Xmax) and a (Xmax) are strongly dependent 
on the value of Xcut- On the other hand, their difference introduced in Eq. (Q above, 
(Xmax) — CT {X^nax), IS rather insensitive to the cutting procedure used. Therefore, the 
possible biasses of the tail of the Xmax distribution do not influence this observable. 
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Figure 7. (Color online) Difference of the observable ycut, evaluated from the Xmax 
distribution truncated at Xcut^ and j/oo, evaluated from the unbiased distribution 
(given by QGSJETOl model for primary protons with energy 10^* eV and 2 • 10"'^^ 
eV). When energy increases, the observed dependence shifts towards the higher Xcut, 
proportionally to the increase of (Xmax) with energy, i.e., the dependence of Ay on 
Xcut — {Xmax) oo remains roughly the same for all primary energies. 



In conclusion, we argue that the spectacular energy dependence of the shower 
maxima distribution reported by the Auger collaboration pTj is not necessarily (or not 
only) due to the changes of chemical composition of primary cosmic rays. The observed 
effect (or, at least, a substantial part of it) seems rather to be caused by the unexpected 
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changes of the depth of first interaction, Xi. However, the energy dependence of {Xi) 
can be affected by a rapid increase of cross section and/or increase of inelasticity in 
energies above 2 ■ 10^^ eV. Both possibihties require an abrupt onset of some "new 
physics" in this energy region and are therefore questionable. The HiRes data [2], where 
the Xmax distribution was truncated and, after that operation, is consistent with the 
proton spectrum, brings in the possible role of biases of the X^ax distribution indicating 
that the ways of analyzing CR data of highest energy still remains an open problem. 
We argue that it would be highly desirable to analyze the observable (Xmax) — o" (Xmax) 
(cf. Eq. (Q) in which fluctuations of the depth of the first interaction, as well as the 
possible biases of the tail of Xmax distribution, are strongly suppressed. This observable 
still depends on the model of multiparticle production and is sensitive to the chemical 
composition of the primary CR. 

Summarizing, though the problem of the chemical composition of CR seem still 
unresolved, we expect that the large spread observed in the results can find an 
explanation by comparing the different sensitivities to the composition of CR in various 
observables. A deeper understanding, both of the hadronic interaction models and of 
the systematics of data analysis, which could bias the results, is also needed. Such a 
detailed analysis, with detailed simulation studies, is, however, outside the scope of this 
paper. 

Acknowledgment 

Partial support (GW) of the Ministry of Science and Higher Education under contract 
DPN/N97/CERN/2009 is gratefully acknowledged. We would like to thank warmly dr 
Eryk Infeld for reading this manuscript. 

References 

[1] Abraham J et al. (Pierre Auger Collaboration) 2010 Phys. Rev. Lett. 104 091101 
[2] Abbasi R U et al., (The High Resolution Fly's Eye Collaboration) 2010 Phys. Rev. Lett. 104 
161101 

[3] Schwarzschild B 2010 Phys. Today 63 15 

[4] Ostapchenko S S 2006 Nucl. Phys. B Proc. Suppl. 151 143 

[5] Alvarez-Muniz J, Engel R, Gaisser T K, Ortiz J A, and Stanev T 2002 Phys. Rev. D 66 033011 

and 2004 Phys. Rev. D 69 103003 
[6] Risse M 2004 Acta Phys. Polon. B 35 1787 

[7] Ulrich R, Bliimer J, Engel R, Schiissler F, and Unger M 2009 New J. Phys. 11 065018 
[8] Wlodarczyk Z 1993 J. Phys. G 19 L133 

[9] Hansen P M, Alvarez-Muniz J, and Vazquez R A, 2011 Astrop. Phys. 34 503 
[10] Ulrich R, Engel R, and Unger M, 2011 Phys. Rev. D 83 054026 
[11] Pirog T, and Werner K 2008 Phys. Rev. Lett. 171101 
[12] Portugal L, and Kodama T 2010 Nucl. Phys. A 837 1 

[13] Shabelski Yu M, Weiner R M, Wilk G, and Wlodarczyk Z 1992 J. Phys. G 18 1281 

[14] Bias de Deus J, Esprito Santo M C, Pimenta M, and Pajares C 2006 Phys. Rev. Lett. 96 162001 

[15] Lindsley J, and Watson A A 1981 Phys. Rev. Lett. 46 459 



On the chemical composition of cosmic rays of highest energy 



12 



[16] Ungcr M (for the Pierre Auger Coll.) 2009 Nucl. Phys. B (Proc. Suppl.) 190, 240 
[17] Kalmykov N N, and Ostapchenko S S 1993 Phys. At. Nucl. 56 346 



